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Summary 

A stochastic design optimization (SDO) methodology has been developed to design components of an 
airframe structure that can be made of metallic and composite materials. The design is obtained as a 
function of the risk level, or reliability, The design method treats uncertainties in load, strength, and 
material properties as distribution functions, which are defined with mean values and standard deviations. 
A design constraint or a failure mode is specified as a function of reliability p. Solution to stochastic 
optimization yields the weight of a structure as a function of reliability p. Optimum weight versus 
reliability p traced out an inverted-S-shaped graph. The center of the inverted-S graph corresponded to 
50 percent (p = 0.5) probability of success. A heavy design with weight approaching infinity could be 
produced for a near-zero rate of failure that corresponds to unity for reliability p {ox p = 1). Weight can be 
reduced to a small value for the most failure -prone design with a reliability that approaches zero (p = 0). 
Reliability can be changed for different components of an airframe structure. For example, the landing 
gear can be designed for a very high reliability, whereas it can be reduced to a small extent for a raked 
wingtip. The SDO capability is obtained by combining three codes: (1) The MSC/Nastran code was the 
deterministic analysis tool, (2) The fast probabilistic integrator, or the FPI module of the NESSUS 
software, was the probabilistic calculator, and (3) NASA Glenn Research Center’s optimization testbed 
CometBoards became the optimizer. The SDO capability requires a finite element structural model, a 
material model, a load model, and a design model. The stochastic optimization concept is illustrated 
considering an academic example and a real-life raked wingtip structure of the Boeing 767-400 extended 
range airliner made of metallic and composite materials. 

Introduction 

Engineers have recognized the existence of uncertainty in material properties, in load, and in 
structural analysis as well as in design constraints. Consider for example the yield strength of a steel that 
is required to design a steel structure. Strength is measured in the laboratory from tests conducted on 
standard coupons. It is commonly observed that repeated tests yield different values for the strength of 
steel. The test data can be processed to obtain a nominal or mean value and a dispersion range or a 
standard deviation. The nominal strength along with a safety factor is used to define allowable strength 
for traditional deterministic design calculations. Alternatively, strength can be considered as a random 
variable with a mean value and a standard deviation. The experimental data can be processed to obtain a 
probability distribution function such as, for example, the commonly used normal distribution function that 
is defined by a mean value and a standard deviation. This concept for strength can be extended to Young’s 
modulus, Poisson’s ratio, density, and so forth, and a probabilistic material model can be generated. The 
procedure can be repeated and a probabilistic load model can be developed for mechanical, thermal, and 
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initial deformation loads. Likewise, a probabilistic design model can be developed for sizing variables like 
depth and thickness of a beam. 

Because of the stochastic nature of load, material properties, and sizing variables, structural response 
consisting of stress, strain, displacement, and frequency become random parameters with mean values and 
standard deviations. The cumulative distribution concept can be utilized to estimate the value of a 
response parameter for a specified level of probability. For example, the value of displacement in a 
particular structure at a location in a direction can be estimated to be less than 0.18 in. for a 25-percent 
probability of success or reliability. The value can change to less than 0.23 in. for 75 percent reliability. 
The concept illustrated for displacement can be extended to failure modes or design constraints, which 
become a function of reliability. In other words, a structure can be designed for a specified reliability 
between 0 and 1 . High reliability can lead to a heavier design. The design is likely to be lighter when 
reliability is compromised. In other words, the weight of a structure becomes a function of the reliability. 

It would be shown that reliability versus weight traced out an inverted-S-shaped graph. 

Reliability-based design optimization requires a probabilistic analysis tool. Several such tools are 
discussed in References 1 and 2. Here, the Fast Probabilistic Integrator (FPI) module of the NESSUS 
code (Ref. 3) is used for probabilistic calculation for the industrial strength wingtip problem. A quadratic 
perturbation method (Ref. 1) is used for the academic problems. The probabilistic response is used to 
formulate the stochastic design problem. It is solved using an optimization testbed, which in the literature 
is referred to as CometBoards (Ref. 4). The probabilistic analysis and design concepts are illustrated for 
an academic example and for an industrial strength aircraft wingtip structure made of metallic and 
composite materials. The subject matter of the paper is presented as follows: Probabilistic analysis is 
covered in the following section, then design optimization is discussed, and conclusions are presented at 
the end. 


Probabilistic Structural Analysis 

Popular probabilistic analysis formulations include Monte Carlo simulation, sampling and stratified 
sampling techniques, the Latin hypercube technique, response surface method, and others. Monte Carlo 
simulation (Ref. 5) is a powerful numerical approach, but it is repetitive and computationally expensive. 
Numerical integration, second-moment analysis, and stochastic finite element methods (Ref. 6) are also 
available. The perturbation method has been used extensively in developing the stochastic finite element 
method because of its simplicity, efficiency, and versatility. Cambou (Ref. 7) proposed a first-order 
perturbation method for the finite element solution of linear static problems. Cornell (Ref. 8) introduced a 
second-moment concept. During the 1980s, the method was developed further by Hisada and Nakagori 
(Ref. 9) for static and eigenvalue problems. The perturbation method was also adopted by Handa and 
Anderson (Ref. 10) for static problems of beam and frame structures. Applications of perturbation 
methods have also been reported by Benaroya (Ref. 11), Elishakoff (Ref. 12), and Schueller (Ref. 13). In 
this paper, for academic problems a quadratic perturbation technique is employed to calculate the mean 
value and the covariance matrix in closed form for stress and displacement. For the wingtip industrial 
problem, fast probability integration (FPI) is used, which is based on a most probable point (MPP) 
concept. The majority of the probability is concentrated near the MPP. The MPP is located and different 
probabilistic methods are used to determine the probability in the failure prone region. The probabilistic 
response is generated from a few deterministic solutions. 

Perturbation Method 

The perturbation method for probabilistic analysis has been developed for both force and 
displacement methods. The method of force used is referred to as the integrated force method (IFM) 

(Ref. 14). The dual of the primal IFM, or the dual integrated force method (IFMD) (Ref. 15), became the 
displacement method. The perturbation technique yields the mean value, and the covariance matrix for the 
response variables like internal force and displacement. The cumulative distribution concept is applied to 
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determine the response for a specified distribution function and a probability level. The parameters for 
perturbation, or the primitive random variables, are separated into three groups. These are loads, material 
properties, and design variables. The load vector encompasses mechanical load, thermal load, and settling 
of support or initial deformation load. It is defined in terms of a mean load vector and a covariance 
matrix. Material properties considered are elastic modulus, Poisson’s ratio, coefficient of thermal 
expansion, and density. Their mean values and the associated covariance matrices are the stochastic 
parameters. The sizing design parameters are the areas of bars, moment of inertia of beams, thickness of 
plates, and so forth. Their stochastic parameters are defined through mean values and a covariance matrix. 
The geometrical parameters like length of a truss bar and spans of a plate are considered as deterministic 
variables. 


Quadratic Perturbation Technique 

In the perturbation technique, the mean value and the covariance matrix of a response variable are 
obtained in two basic steps. First, a response variable is expanded in a Taylor’s series with respect to the 
primitive random variables. The linear and quadratic terms are retained in the series expansion. An 
application of the expectation operator E[..] yielded the expressions for mean value and the covariance 
matrix. 

Prior to the expansion, the primitive parameters are scaled to obtain a normalized variable (« q ) with a 
zero mean and a small variance. Consider, for example, the coefficient of thermal expansion a with a 
mean \i a and a standard deviation a a . Its normalized variable q a is defined as 


4a 


a ~Pq 

Pa 


(la) 


or 


a = M 1 + ?a) 

The mean value of q a is obtained by using the expectation operator (E[..]): 
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The variance of a is assumed to be much smaller than the square of its mean. This justifies the use of 
a Taylor expansion of response variables in terms of the normalized primitive random variables. The 
Taylor series expansion is discussed considering the example of the internal force vector {F}. The 
expansion can be expressed as 

{F} = {F} + £{F} j + ^ {#} j + higher order term (3) 


where 


and 


(faf £(•)) = &• 


5 (.) 




'{?}={«} 


({qf C(')C T {?}) = £?/?,• 

Uj 


FIT 

dqfiqj 




Similar expansions are also obtained for the governing matrix [S] of IFM, the load vector {F}, 
displacement { X }, the flexibility matrix [G], and initial deformation {(3}, where their deterministic (or 
over-bar) quantities are calculated by setting the normalized primitive random variables to zero (q t = 0). 
Four steps are followed to derive the expressions for the mean value and the covariance matrix of the 
force vector: 

Step 1. Substitute the Taylor series expansions into the IFM governing equation, see 
References 1 and 2. The left-hand side contains a single zero-order term, two linear terms, and three 
quadratic terms in the normalized primitive random variables, and the right-hand side has one term of 
each order. 

Step 2. Equate equal order terms. The zero-order term gives the deterministic solution. The linear 
terms drop out. The second-order terms produce an equation with four quadratic terms on the right-hand 
side and a single term on the left-hand side. 

Step 3. Take the expectation E[..] of terms in Step 2 to obtain the mean value of the force vector. 

Step 4. Take the expectation E[..] of the product {F} {F} r , and subtract the square of the mean of the 
force vector, as calculated in Step 3, to obtain the covariance matrix of the force vector. 

The expressions for mean value and covariance matrix for force and displacement vectors via IFM, 
listed in Reference 1, is not repeated here. The process was repeated next for the IFMD to obtain similar 
expressions, also listed in Reference 1. The derivation included mechanical and initial deformation loads. 
The mean value of the force and displacement contained zero-order and quadratic terms, but not linear 
terms, in the Taylor series. However, the covariance matrix retained the zero-order, linear, and quadratic 
terms of expansions. For simple academic structures the expressions were programmed in closed form in 
Macsyma software (Ref. 16). 
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Illustrative Example: Three-Bar Truss 


Stochastic response analysis is illustrated considering the example of a three-bar truss shown in 
Figure 1. The geometrical dimensions of the structure like the coordinates of the four nodes and bar 
lengths are considered deterministic in nature. The structure has a total of 10 primitive random variables 
for the material properties, load, and sizing design parameters. 

Material Model: The truss is made of steel, and it has two random material variables: the Young’s 
modulus and the coefficient of thermal expansion. The Young’s modulus E has a mean value of 
\i E = 30 000 ksi and standard deviation of 3000 ksi. The coefficient of thermal expansion a has a mean 
value of p a = 6.6 x 10“ 6 /°F and standard deviation of 6.6x1 0“ 7 /°F. 

Mechanical Load Model: The structure is subjected to mechanical and thermal loads as well as loads 
due to settling of support. Mechanical load is defined by its mean value and covariance matrix. The load 

| is applied at the free node 1 . 



Mean value: 


w 



Covariance matrix: cov 


'\Px\ 

' [25.00 

6.25 “ 

[\ p y J 

, 6.25 

25.00 


Thermal Load Model: Node 1 of the structure is subjected to a change of temperature, with a mean 
of 1117 = 50 °F and a standard deviation of 5 °F. The other three nodes are at ambient temperature. Thermal 
load is obtained for the bar temperature, which is calculated from the average of the temperatures at its 
two connecting nodes. 

Load Due to Settling of Support: Support node 2 settles in both the x- and y-directions by amount 
| inches as shown in Figure 1. Settling is accommodated by a mean value and a covariance matrix. 




Figure 1 . — Three-bar truss, where P is load, F is force, and 
A is setting of support. 
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Mean value: 



in. 


Covariance matrix: cov 


r 

^2X 

A 


'25.0 37.5 ' 

v 1 

[^27 J 

^ ) 


37.5 225. 0_ 


xlO“ 


Design Model: The areas of the three bars (A u A 2 , and^4 3 ) are considered to be the design variables. 
Their mean value and covariance matrix are as follows: 
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Covariance matrix: cov 
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The 10 random variables of the problem are scaled to obtain the 10 primitive variables q i, qi,...,q io, 
which are defined as 


A 

= ^ 1 ( 1 +^ 1 ) 

- Oi (l + <h ) 

A 

= M 1+ *J 

= \y 2 (\ + q 2 ) 

A 

= ^ 3 ( 1+ ^ 3 ) 
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E 

= (l + Qe) 
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a 

= M 1 + ?«) 
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Px 


= 06( 1 + ?6) 

Py 


= M 1 + ? 7 ) 
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A 

= ^A, (l + « Al ) 

= H 9 ( 1 + ^ 9 ) 

a 2 

= ^a 2 (l + ? A2 ) 

= 0io( 1 + #io) 


The normalized primitive random variable, {q}, with zero mean, and standard deviation given by the 
ratio of the standard deviation to the mean of the corresponding stochastic variable is assumed to be of 
order o(l). This justifies the use of a Taylor series expansion in {q}. Only up to quadratic terms are 
retained in the series. 
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The deterministic response for forces and displacements are as follows: 


Bar forces: {F} = {p^ } = 


62.76 

61.24 

-7.95 


kip 


r . ( , f 0.201 ] 

Nodal displacements: \X) = |p x j = \ q [ * n ’ 


The probabilistic response for bar forces consists of the mean values and the covariance matrix as 
follows: 


Mean values of bar forces: jp^} = 


62.78 

61.22 

-7.93 


kip 


Covariance matrix of forces: cov(F) 


19.27 

1.78 

-4.48 

1.78 

18.93 

-8.83 

-4.48 

-8.83 

21.78 



n 2 


1.0 

0.093 

-0.219' 

Correlation matrix for forces: 


“P F~ 

0.093 

1.0 

-0.435 


_(Pi)(Py)J 

-0.219 

-0.435 

1.0 


Likewise, the probabilistic response for nodal displacements are 


Mean values of displacements: jp^} 


0.197 1 

\ in. 
— 0.237 J 


Covariance matrix of displacements: cov(x) 


12.2 

-8.9 


-8.9 

8.7 


xlO -4 


Correlation matrix for displacements: p x = 


1.0 

- 0.866 


- 0.866 

1.0 


where the elements of and p^are the correlation coefficients. 

For the bar forces, there is very little difference between the mean value {p^j and the deterministic 
solution {p F } . For the nodal displacements also, there was little difference between the mean value {p^j 

and the deterministic solution {p x } . For all examples, big or small, it was observed that the mean value was 

quite close to the deterministic solution, which however, has not been proved analytically. 

The standard deviations for the bar forces were about 7 percent for bars 1 and 2. However, it was 
59 percent for the bar 3. The variation in the primitive variables like mechanical, thermal, and settling of 
support loads was between 5 and 10 percent. The variation in the response or bar forces of the basic 
structure is comparable to that for the load variation. However, the variation in the redundant bar force at 
59 percent is 6 times that of the load. 
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Figure 2. — Probability density function for bar 1 force and vertical load. 



Figure 3. — Cumulative distribution function for bar 1 force 
and vertical load. 


The standard deviations for the displacements are between 12 and 17 percent, which are higher than 
the load variation. The variations for the displacements are more pronounced than that for the bar forces 
because these are susceptible to loads as well as to material properties and design variables. 

The probability density function and the cumulative distribution function for the bar 1 force and 
magnitude of the vertical load are depicted in Figures 2 and 3, respectively. There is 50-percent likelihood 
that the first bar force (with a mean of 62.67 kips) is between 59.80 and 65.73 kips. 

The value of the response variables for percent probability of success p = 50, 25, and 75 percent are 
listed in Table I. A 5-percent change was seen in the forces of bars 1 and 2 between the p = 50 percent and 
p = 15 percent levels. For the bar 3 force, the change noted was 40 percent. A 12-percent change was seen in 
the horizontal displacement between the/) = 50 percent and p = 15 percent levels. It was only 8 percent for 
the vertical displacement. 


TABLE I.— STOCHASTIC RESPONSE FOR THREE-BAR TRUSS 


Response variable 

Upper limit of response variable (range) for 
probability of occurrence,/? 

50 percent (mean) 

25 percent 

75 percent 

Bar force, kips 




Fx 

62.76 

59.80 (95%) 

65.73 (105%) 

F 2 

61.24 

58.30 (95%) 

64.17(105%) 

F 3 

-7.95 

-4.80 (60%) 

-11.09(140%) 

Displacement, in. 




u 

0.201 

0.178 (88%) 

0.225 (112%) 

V 

-0.241 

-0.221 (92%) 

-0.261 (108%) 
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The sensitivity of the bar 1 force and the vertical displacement at the 75 -percent probability of 
occurrence level with respect to the 10 primitive random variables are shown in Figures 4 and 5, 
respectively. The maximum value along the y-axis in Figures 4 and 5 are normalized to unity. The bar force 
is equally sensitive to the areas of bars 1 and 2 in magnitude, but opposite in directions. The sensitivity of 
the bar 1 force with respect to the bar 3 area is very small. The vertical displacement is equally sensitive to 
the areas of bars 1 and 2, but is insensitive to the bar 3 area. The Young’s modulus has little effect on the 
bar force but has the most effect of all the primitive variables on the vertical displacement. The 
sensitivities with respect to thermal coefficient and temperature are identical because the analysis uses 
their product as a single entity. Both load components strongly effect the bar force, whereas only the 
vertical load effects the vertical displacement to the same degree. Temperature and settling of support 
loads have a moderate effect. 



Figure 4. — Sensitivities for bar 1 force at 75 percent probability level with respect to 
10 random variables. 



Figure 5. — Sensitivities for vertical displacement at 75 percent probability level with 
respect to 10 random variables. 
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Comparison With Other Methods 

Stochastic response via perturbation method was compared with other probabilistic formulations 
considering the example of the three-bar truss shown in Figure 1 with the same 10 random variables for 
the material properties, load, and sizing design parameters. A summary of response calculated by the four 
methods is depicted in Table II. Two different computers were used in the calculations shown in Table II. 
A Dell Inspiron desktop with four central processing units (CPUs) and 3.2 GHz processor was used to 
obtain most of the results. The FPI code shown in the last column was run in a very fast Red Hat Linux 
dual CPU workstation with a 32 GHz processor. The four methods compared were 

(1) Perturbation method (PM) 

(2) Direct Monte Carlo simulation (DMCS) 

(3) Latin hypercube simulation (LHS) 

(4) Fast probabilistic integration of the NESSUS code (FPI) 


TABLE II.— PROBABILISTIC RESPONSE COMPARED FOR THREE-BAR TRUSS 


Parameter 

Perturbation method 
(PM) 

Direct Monte Carlo 
simulation (DMCS) 

Latin hypercube 
simulation 
(LHS) 

Fast probability integrator 
(FPI) 

Mean 

value 

Standard 

deviation 

Mean 

value 

Standard 

deviation 

Mean 

value 

Standard 

deviation 

Mean 

value 

Standard 

deviation 

Force : < 

F 2 
A 

kip 

> 


62.771 

61.24 

-7.95J 



4.39 

4.35 

4.67 



f 62.77 
61.27 

[-7-93. 

r* 


4.39' 

4.35 

4.67 



[62.761 

61.25 

[-7*96 J 



4.39 

4.35 

4.67 


< 

62.78 

61.21 

-7.93 



4.39 

4.35 

4.67 


Stress: < 

g 2 

ksi 

> 


63.291 

61.70 

-3. 98 J 

> 


6.71 

6.20 

2.34 



f 63. 1 6 
61.58 
[-3.97^ 

> 


5.89 
5.42 
2.7 9 



[63.151 

61.57 

[-3.99J 



5.77 

5.38 

2.75^ 



62.78" 

61.21 

-3.96 

> 


6.34' 

6.03 

2.69 


M in ' 

Displacement : < S- 

w 

- 

\ 0.20 
l” 0 - 24 . 


i 

0.004 

0.003 

} 


[ 0.20] 
[-0.24] 

1 

1 

0.004 

0.003 

} 


[ 0.20] 

[-0.24] 

1 

i 

0.004 

0.003 

} 

< 

f 0.20' 

l-°- 24 


I 

0.004 

0.003 

} 

Computation time | 

CPU seconds 

7.0 

4245 

390 

1 

Normalized time 

1.0 

606 

56 

___ 


The stiffness method as implemented in the ANSYS code (Ref. 17) was used in the DMCS and in the 
LHS. Response was calculated via the primal and the dual integrated force methods (IFM and IFMD, 
respectively) in the perturbation and fast probability integration techniques. Both IFM and IFMD yield 
identical solutions for deterministic as well as stochastic calculations even though the structure of 
equations differed, at least in appearance (Ref. 1). Response for the three bar forces, stresses, and 
displacements by the four methods (PM, DCMS, LHS, and FPI) and time to solution (CPU seconds) are 
depicted in Table II. For bar forces, the mean values and standard deviations were in good agreement for 
all four methods. The displacements showed an almost perfect match across the four methods. There was 
a minor mismatch among the methods for bar stress. The maximum difference was about 0.8 percent in 
the mean value of the stress, whereas it was about 19 percent for standard deviation for the redundant 
member. Direct Monte Carlo simulation required 12 500 samples for convergence, whereas 1000 samples 
were sufficient for the Latin hypercube simulation. The time to calculate the response was very small, 
between 1 to 7 CPU seconds by the perturbation method as well as by the fast probability integrator. The 
calculation time increased many times for the Monte Carlo and Latin hypercube simulations. Monte Carlo 
simulation required about 4245 s, which corresponded to 606 times that required by the perturbation 
method. Latin hypercube method took 56 times as long. Overall, the performance was satisfactory for all 
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four methods. The Ph.D. dissertation of Wei (Ref. 2) provides merits and limitations of different 
probabilistic solution methods. 


Probabilistic Design Optimization 

Probabilistic design optimization accounts for the uncertainties in the model, like those in design 
variables, in loads, in strength, and in material properties. The design is obtained as a function of the risk 
level, or reliability,/?. A design with a 50-percent rate of success (corresponding to one failure in two 
samples, or A = 2 and p = 0.5) is the mean- valued design. For the most reliable design the parameter p 
can approach unity. For example, for one failure in 50 000 samples,/? = 1 - 1 /N= 0.99998. For a design 
that is not reliable or prone to failure,/? can approach 0, (e.g., for 999 failures in 1000 samples, 
p = 0.001). The constraints that represent failure modes are defined for a specified rate of failure. 

A probabilistic design optimization problem can be formulated with two merit functions: 

(1) Minimize the mean value of the weight 

(2) Minimize the variation or the standard deviation of the weight 

Consider for example, the three-bar truss shown in Figure 1. The area or diameter of a bar becomes 
the design variable. Probabilistic calculation can consider the determination of mean value and a standard 
deviation in the diameter. However, the standard deviation in the diameter is a manufacturing issue. A 
design process uses available members with specified diameter or cross-sectional areas. 

Probabilistic optimization can be separated into reliability-based design and robust design. 

Reliability -based design minimizes the mean value of the objective function. Robust design 
simultaneously minimizes the mean value as well as the standard deviation of the objective function. This 
paper emphasizes reliability-based optimization. However, for completeness, the robust design concept is 
discussed for the three-bar truss shown in Figure 1. 

The mean value and standard deviations of the objective function (f) can be defined as 

Mean value: = E[f (x)] = JJ...Jf(x)p(x 1 ,X 2 ,...x^)dx 1 dx 2 ...dx^ (5a) 


Standard deviation: Gy = E 


2 2 

(f(x)-p/) =JJ...J(f(x)-p / ) p(x 1 ,x 2 ,...x /t )dx 1 dx 2 ...dx /i: (5b) 


where 

p(xi, x 2 ,. . .Xk) is the joint probability density function 
f(x) is the merit (or weight) function 
X; are the design variables 
E[..] is the expectation operator 

It is difficult to define the joint probability density function for an industrial-strength structural design 
problem. Here, normal distribution is assumed for the random variables, which are further considered to 
be independent. The robust design optimization problem can be reformulated as follows: 


Minimize both mean value py and standard deviation ay / = [^py,GyJ 
E[g z - (x)] + (x)) < 0 within prescribed lower and upper bounds x L 


for constraints 

< X < Xjj 
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where 

f is a vector of two objective functions 

p/is the mean value 

<jf is the standard deviation 

P z is a constraints (g z ) satisfactory index 

Pareto optimality criteria (Ref 18) is a classical technique to solve a multiobjective problem. Because 
of difficulty, however, the problem is transformed to a single composite objective function and solved. A 
composite objective is obtained by a linear combination of mean value and standard deviation with a 
weighting factor a in the range 0 < a < 1 . The function of merit F can be written as 

F = (l-„)!t + „i (6) 

!% G w 

* * 

where \x w and <5 W are the mean and standard deviation of weight, respectively, when solved 
independently. 

The robust design of the three-bar truss shown in Figure 1 uses the data given for its analysis, along 
with a weight density that has a mean value of 0.289 lbf/in. 3 with a standard deviation of 0.005 lbf/in. 3 . 
The mean value and the standard deviation for the strength are 20 and 2 ksi, respectively. The optimum 
solution is given in Table III for different values of the parameter a for a 50-percent probability (p = 0.5). 
Mean value and standard deviation versus the parameter a are plotted in Figure 6 with an exaggerated 
scale along the y-axis. 


TABLE III.— ROBUST DESIGN FOR THE THREE-BAR TRUSS 


Design parameters 

Weight factor, 
a 

0 

0.25 

0.50 

0.75 

1.00 

Area, in. 2 : ■ 

^2 

kJ 



1.1042 

0.7359 

0.1000 



1.1042 

0.7359 

0.1000 



1.1042 

0.7360 

0.100 



1.1042 

0.7360 

0.1000 



1.1042 

0.7361 

0.1000 


f mean ] 
Weight, lbf: \ > 

[ variance J 

1 

70.4908 

8.4504 


1 

70.4889 

8.4502 


i 

70.4870 

8.4500 


1 

70.4863 

8.4499 


1 

70.4870 

8.4500 



The mean value and standard deviation versus the parameter a traced out the typical graph as shown 
in Figure 6. The robust design corresponding to the intersection point MS represents a simultaneous 
minimum for the mean value of weight as well as its standard deviation. The intersection occurred for 
a = 0.58. Weight has a mean value of 70.488 lbf and a standard deviation of 8.45 lbf. Critics may label 
the robust design optimization calculations as academic because the range of variation is very small. The 
variation in the mean value of the weight was in the range 70.4863 to 70.4908 lbf. The range for standard 
deviation was 8.4499 to 8.4504 lbf. Robust design optimization is discussed further in Reference 2. 
Reliability-based optimization, also referred to as stochastic design optimization, is discussed next. 
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70.492 



70.484 -I 1 1 1 1 

0.00 0.25 0.50 0.75 1.00 

Weighting factor, a 

Figure 6. — Robust design optimization for three-bar truss. 
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Stochastic Design Optimization (SDO) 

A stochastic optimization methodology has been developed to design components of an airframe 
structure that can be made of metallic and composite materials. The design is obtained as a function of the 
risk level, or reliability,/). The design method treats uncertainties in load, strength, and material properties 
as distribution functions, which are defined with mean values and standard deviations. A design constraint 
or a failure mode is specified as a function of p. Solution to stochastic optimization yields the weight of a 
structure as a function of p. Optimum weight versus p traced an inverted-S-shaped graph. The center of 
the S graph corresponds to 50 percent,/? = 0.5, probability of success. A heavy design with weight 
approaching infinity can be produced for a near zero rate of failure that corresponds to unity for p, or 
p= 1. Weight can reduce to a small value for the most failure-prone design (p = 0). Reliability can be 
changed for different components of an airframe. For example, the landing gear can be designed for very 
high reliability, and reliability can be reduced to an extent for a raked wingtip. The SDO capability is 
obtained by combining three codes: (1) The MSC/Nastran code is the deterministic analysis tool, (2) the 
fast probabilistic integrator, or the FPI module of the NESSUS software, is the probabilistic calculator, 
and (3) NASA Glenn Research Center’s optimization testbed CometBoards became the optimizer. The 
SDO capability requires four models: a finite element structural model, material model, load model, and 
design model. This paper illustrates the design capability by considering the example of a Boeing 
767-400ER raked wingtip. 


Formulation of Reliability-Based Design 

The design testbed CometBoards has been extended into the stochastic domain. The objective is to 
determine the mean values of the design parameters that optimize the mean value of a merit function 
(such as weight), subject to a set of constraints. Variance in the weight can be back calculated. For 
stochastic optimization, the constraints gj are derived from random response variables within prescribed 

random upper and lower bounds g^ and g L j , respectively, with a specified percent probability 

(p between 0 and 1), as: P(gj < gj < g^)>p. A stochastic constraint (e.g., a stress limitation x < x 0 ) can 
be expressed as 
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Standard deviation 


Stochastic 

g(p) 


Mean value 

i- 



-l 

+ 


^ 10 




<p» 


I 


2 2 

<3 r +G2. 

T 1 T 10 


<0 


( 7 ) 


where 

p is the probability of occurrence 
p Ti is the mean value of stress 

p Tio is the mean value of strength (stress allowable) 

cp* is the parameter phi critical 

p Ti is the standard deviation of stress 

cj t is the standard deviation of stress allowable 

x io 


The first term with superscript “Mean value” resembles a deterministic stress constraint. It is written 
in terms of the mean values of the response variable, here stress Xi has a mean value p Ti . Strength or 

allowable stress Xi 0 has a mean value of \i T ^ . The second term with superscript “Standard deviation” 

accounts for the standard deviations a Ti and a Tio in stress and strength xi and Xi 0 , respectively. The 

inverse of the cumulative distribution function for the standard normal at probability level p is cp *(/?). The 
stochastic constraint reduces to the familiar deterministic form, when the variation is set to zero. 



a T ^0 

T 1 



Lim< 

a T 0 




T 10 

u T —> allowable 

1*40 

\ 

.^0 


-1 


+ <P*(p) 




2 2 

<s T +a;r 
T 1 T 10 


<0 


Deterministic constraint 


f 

\ 


1 


— 1 

{^o 

J 


<o 


( 8 ) 


The derivation of the constraint expression in Equation (7) is straight forward. Consider a stress 
constraint (|x| < r 0 ) with a probability p of occurrence of p: 


p(\x\ZT 0 )zp 

p(|t|-t 0 <0)>/7 


Considering the case when x is greater than zero, the difference between the two random variables is 
set to a new random variable iR: 


iR = x-x 0 (10) 

The new random variable is normalized to obtain a standard normal random variable O with mean 0 and 
variance 1: 


0 = 






( 11 ) 


Thus 
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< o) > p 

V G « CT 9? ) 


The probability p(0 < x) is the definition of the cumulative distribution function for the standard 
normal F 0 . Thus 


F® >P 

v 


The minimum value of O at which the probability level, p , is satisfied is obtained from the inverse of 
the cumulative distribution function for the standard normal. This value is labeled “cp 

<P*(p) = F®(/>) 

*i \ m? (14) 

cp (p)<-^ 




c 9? - c i + G x 0 


cp*(p)<- 


m-i\, 


( P*( J P)V C ? +G ? 0 ^-^T+^Xo 

^x-^x„ + 9 *(P)\[<^+°L 


Finally, taking the case x < 0, normalizing with respect to the mean of the allowable and combining 
gives 


^ ll+ 

/ 2 2 
* < \ r* + 

9 yp) 

_ ^0 J 

* ^0 


or g < 0 
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As mentioned earlier, the constraint has two parts. The first part is similar to a deterministic constraint 
specified in terms of the mean value of the variable. The second part contains the contribution from 
stochastic that included the inverse function and square of standard deviations as well as the mean value. 
For p = 0.50, cp is zero, and the constraint degenerates into a standard deterministic constraint. Observe 
that jli x is a response variable that is obtained from stochastic calculations and is not the same as the 
solution to the deterministic system. 

In the above derivation, the two cases x > 0 and x < 0 were treated separately. In doing so, an 
assumption was made that p T » g t . This is reasonable when the constraint is active. For p = 0.25, cp* is 
-0.6745, and the constraint is relaxed. If it were an active constraint, the stochastic process becomes 
equivalent to increasing the allowable value. This would have a tendency to generate a lighter optimum 
design. For p = 0.75, cp* is 0.6745, the constraint is tightened. If it were an active constraint, the stochastic 
process becomes equivalent to decreasing the allowable value. This would have a tendency to generate a 
heavier optimum design. The stochastic design will have a tendency to produce a lighter design when p is 
less than 50 percent, but a heavier design otherwise. The formulation of reliability-based design 
optimization is quite similar to that for the deterministic problem. The mean value of the design variable 
is generated to minimize the mean value of weight subjected to the stochastic constraints defined in 
Equation (7). 


Numerical Examples for Reliability-Based Design Optimization 

Reliability -based optimization is illustrated considering two components of an airframe structure. The 
first component is a web of an airframe stabilizer, and the second component is a racked wingtip. 
Examples are proprietary components of the Boeing Company and provided to us as a courtesy to 
advance reliability-based optimization concept for industrial applications. Sufficient information will be 
given to illustrate the concept. 

Design Example 2: Web Plate of an Airframe Stabilizer 

The inverted-S-shaped graph in reliability-based design optimization is illustrated considering the 
example of a web plate of an airframe stabilizer. The steel web plate, which was part of an aircraft 
stabilizer structure, was about 100 in. long, 12 in. wide, and 0.072 in. thick. It was idealized as a planar 
structure in the x,y-plane. The finite element model had CQUAD4 and CTRIA3 elements. Load was 
applied in the x,y-plane. Load was obtained by projecting the reactive load from the stabilizer into the 
plane of the web. The maximum von Mises stress for the initial design was 100 units psi with a maximum 
displacement of A milli-inch. Load, material properties, and design parameters were considered as 
random variables. The mean values were taken equal to their deterministic values with a 10-percent 
standard deviation. For example, the mean value of Young’s modulus was 30 000 ksi, (which was its 
deterministic value) with a 10-percent standard deviation of 3000 ksi. The stochastic response calculation 
used three different probabilistic analysis methods. All three methods used MSC/Nastran for deterministic 
calculation. The FPI module was also used by all three methods. The methods used were 

Method 1: This method used the first probabilistic integrator with normal distribution for all 
variables. 

Method 2: A neural network and a regression method were trained for the three design variables. 
Normal distribution was used with the neural network approximation. 

Method 3: This method replaced normal distribution with the Weibull function in method 2. The 
regression method was used for approximation. 
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There was little difference between the mean value and deterministic value for stress by all three 
methods. The standard deviation was about 9 percent, which compared well with the 10-percent deviation 
in the primitive random variables. The design model had three random variables, consisting of one 
thickness along the web boundary, one in the plate central zone, and one in between (the intermediate 
region). To begin optimization calculations, the mean value of the variables were taken equal to the 
deterministic optimum solutions. The standard deviations were 10 percent of the mean value for all three 
of the design variables. CometBoards testbed was used for design optimization. Optimum weights for 
different probability of failure are given in Table IV. The rate ranged from a vulnerable design (with 999 
failures in 1000 samples) to a reliable design with 1 failure in 1000 samples. 

Weight versus probability of success is plotted in Figure 7. There is some deviation in the weight 
calculated by the three different analysis methods as shown in the figure. For a 50-percent probability of 
success all three methods converged to 3.22 lbf for the optimum weight. For a reliable design with 1 
failure in 1000 samples the first method produced about a 10-percent higher weight than that obtained by 
the neural network technique. The difference reduced to 5 percent for the regression method. 


TABLE IV.— OPTIMUM WEIGHT VERSUS PROBABILITY LEVEL FOR WEB PLATE 


w 

Probability 


Optimum weight of web plate, 


level, 


lb 


(X failures in 

P 

Method 1 

Method 2 

Method 3 

1000 samples) 


FPI and normal 

Neural network and 

Regression method and 



distribution 

normal distribution 

Weibull function 

(999) 

0.001 

1.16 

1.78 

1.75 

(975) 

0.025 

1.89 

2.26 

2.22 

(900) 

0.1 

2.33 

2.56 

2.54 

(700) 

0.3 

2.84 

2.94 

2.93 

(500) 

0.5 

3.22 

3.22 

3.22 

(300) 

0.7 

3.61 

3.52 

3.52 

(100) 

0.9 

4.21 

3.99 

3.98 

(25) 

0.975 

4.81 

4.47 

4.41 

(1) 

0.999 

5.96 

5.40 

5.16 



Probability level 


Figure 7. — Inverted-S graph for optimum weight versus probability level for 
web plate. 
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The Inverted-S-Shaped Graph 

The weight versus probability of success for the web plate given in Table IV is plotted in Figure 7. 
Weight increased when reliability exceeded 50 percent. Weight decreased when the reliability was 
compromised. The weight versus reliability traced out the inverted-S-shaped graph. A design can be 
selected that depends on the level of risk acceptable to the situation. The inverted-S-shaped graphs were 
also generated for 20 examples given in Reference 2. The illustration in Figure 8 provides a simple 
explanation for the inverted-S shape of the graph. 

Consider a deterministic design calculated for an allowable stress limit of 25 ksi with a 100 lbf weight 
for the structure (see Fig. 8(a)). The structure is redesigned next for a 50-percent probability of success. 
Assuming a 1.5 safety factor in applied loads, the stress can be approximated to 17 ksi and the 
proportioned weight can be approximated as 67 lbf (see Fig. 8(b)). Let us consider a reliable design with 1 
failure in 100 000 samples. The stress value is likely to increase to 24 ksi because of an increase in the 
corresponding area under the distribution function (see Fig. 8(c)). The weight has to be increased to about 
80 lbf because of a high value for stress, as shown in Figure 8(c). Consider next a failure -prone design 
with 90 failures in 100 samples. The stress can be less, like 7 ksi, because of a reduced area under the 
distribution function (see Fig. 8(d)). The weight can be reduced to about 28 lbf because of the low stress 
level. In reliability-based design optimization, weight can become very heavy when reliability approaches 
unity; likewise, a lightweight design can be obtained when reliability is compromised. 


Assume 

stress a = 25 ksi 
weight = 100 lbf 


(a) Deterministic design 



Mean value ► Stress 

(b) Probabilistic design for p = 0.5 


Weight = 80 lbf 1 failure in 



Weight = 28 lbf 


(d) o=7 ksi 


90 failures in 
100 samples 


Stress 


Figure 8. — Basic reliability design concept; a qualitative illustration for induced stress, (a) Deterministic 
design, (b) Design for mean value of reliability p = 0.5. (c) Design for more than mean value, (d) Design 
for less than mean value. 
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Figure 9. — Model of raked wingtip of 
Boeing 767-400 extended range airliner. 


Design Example 3: Raked Wingtip Structure 

Reliability-based optimization is illustrated considering the example of a composite airframe 
component of the Boeing 767-400 extended range airliner, the raked wingtip structure shown in 
Figure 9. The structure is fabricated out of component parts made of metallic and composite materials. 
The available industrial design is referred to as the initial design. It is required to generate a deterministic 
optimum design as well as a reliability -based design for the component. The problem is treated in five 
separate sections: 

(1) Deterministic analysis 

(2) Deterministic optimization 

(3) Probabilistic analysis 

(4) Probabilistic analysis for strain and displacement 

(5) Reliability-based design optimization 

Deterministic Analysis 

The component has a wing-box type of construction with face sheets and web panels made of 
composite materials with an aluminum honeycomb core. Its components such as the root rib, tip rib, and 
leading edge parts are made of aluminum, and a set of titanium members are also used. The MSC/Nastran 
code was used as the deterministic analysis tool. The finite element model had about 3000 nodes and 
17 000 degrees of freedom. The model was made of 3500 elements, consisting of CBAR, CBUSH, 
CELAS1, CHEXA, CONROD, CPENTA, CQUAD4, CROD, C SHEAR, and CTRIA3 elements. For 
aluminum, titanium, composite, and honeycomb materials, typical values were used for the material 
properties like Young’s modulus, Poison’s ratio, shear modulus, and others. The structure was subjected 
to eight load cases. The CPU time for a single MSC/Nastran run was about 5 s for static analysis, and it 
required almost 1 min to calculate a set of 20 frequencies on a Red Hat Enterprise Linux workstation, 
Release 4. Overall response of the initial model for the eight load cases is presented in Table V. 
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TABLE V.— NORMALIZED RESPONSE OF COMPOSITE RAKED WINGTIP 
(MAXIMUM VALUE IS SET TO 100 UNITS) 


Load case 

Strain energy, 
U, 

in. -kip 

Maximum 

displacement, 

in. 

Maximum von Mises 
stress for metal, 
ksi 

Maximum principal 
micro strain in 
composite laminates 

1 

17.22 

-41.15 

43.56 

32.49 

2 

73.70 

84.50 

76.85 

72.96 

3 

63.70 

77.41 

74.91 

65.86 

4 

61.11 

75.73 

75.37 

66.13 

5 

16.48 

-40.50 

40.48 

34.81 

6 

100.00 

100.00 

100.00 

100.00 

7 

82.04 

91.01 

96.69 

93.44 

8 

52.22 

77.49 

70.24 

55.96 


The strain energy U in the first column in Table V was obtained as one-half of the product of load P 
and displacement X over all nodes using the following formula: 


U= — I | (strain x stress) dr 


Volume 


( All nodes and directions 'N 

i 

7=1 


PjXj 


(18) 


The strain energies for the eight load cases were in the range 16 to 100 normalized units. Load case 6 
contained the peak value for strain energy of 100 units. The load case also contained the maximum 
displacement, maximum strain in the composite components, and maximum stress in metals. The 
structure was designed for the critical load case 6 and checked against other load cases. The concept to 
design for a load case with maximum strain energy drastically reduced the number of calculations in 
design optimization, and it was subsequently proved correct for the problem. 

Deterministic Optimization 

The objective of the optimization problem was to reduce the weight of the wingtip without changing 
its geometrical configuration. Only the thicknesses of the components were allowed to be adjusted within 
prescribed bounds. For design optimization the raked wingtip was separated into several substructures 
consisting of metallic and composite components. From a preliminary design optimization study it was 
concluded that the thickness of titanium components and aluminum parts like, root rib center section, 
leading and trailing edge webs, as well as leading edge and tip rib skins cannot be changed from a 
manufacturing consideration. Design variables associated with these parts were considered passive. The 
remaining four composite members were grouped to obtain a total of 13 active design variables. Grouping 
was on the basis of the number of laminates. For example, the first design variable represented a 
proportional thickness parameter for the upper and lower skin panels. The panel was made of laminates 
over a honeycomb core. This variable contained a total of 21 1 CQUAD4 and CTRIA3 elements. 

Likewise, the design variable 10 represented the area parameter for spars in the front web. The areas of 
the spars were grouped to obtain a design variable and 92 CROD elements were used in the discretization. 
Other design variables were defined in a similar manner. 

For constraint formulation, the structure was separated into 203 groups of elements to obtain a total of 
203 strain constraints for the upper and lower skin panels as well as the spars. The composite rod 
elements were grouped to obtain 16 more strain constraints. Three translations and one rotation were 
constrained at a tip node of the structure for load case 6 as well as for load case 7 to obtain eight 
displacement constraints. The design model has a total of 227 behavior constraints. Constraint can be 
imposed on principal strain, for a failure theory, or on a strain component. Displacement limitations along 
x-, y-, and z-directions were set at A 2 , and A 3 in., respectively. Maximum rotation was not to exceed 
0°. Design optimization was performed using the testbed CometBoards. Deterministic optimum solution 
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is given in Table VI. Only a normalized optimum solution can be given because of the proprietary nature 
of the data. 


TABLE VI.— NORMALIZED DETERMINISTIC OPTIMUM SOLUTION 


Design variable 

Weight 

Original design 

Design model 1 

Design model 2 

100 units a 

84.0 units 

80.0 units 

Change in percent 
(100 represents no change) 

1 

100.0 

70.1 

70.1 

2 

100.0 

104.5 

99.83 

3 

100.0 

110.3 

111.83 

4 

100.0 

75.64 

76.97 

5 

100.0 

44.85 

44.45 

6 

100.0 

123.35 

80.96 

7 

100.0 

84.17 

94.45 

8 

100.0 

44.08 

44.88 

9 

100.0 

86.12 

83.59 

10 

100.0 

88.09 

88.09 

11 

100.0 

73.83 

73.83 

12 

100.0 

100.00 

100.00 

13 

100.0 

100.00 

100.00 

Active strain 
constraints 

100.0 

6 active and 
9 nearly active 

6 active and 
9 nearly active 

Displacement in 
z-direction, in. 
(limit was A in.) 

100.0 

Nearly active 

Nearly active 

Rotation, deg 
(limit was 0°) 

100.0 

4.48 

4.72 

Frequency, Hz 

100.0 

16.45 

20.33 


formalized to 100 units. 


Stress and strain obtained from the finite element model were presumed not to be accurate at some 
metal and composite interface locations. Such localized regions were avoided in design calculations by 
excluding the associated strain constraints. Table VI, column 2, model 1, refers to a configuration that 
was obtained by excluding the strain constraints for a set of interface elements. Likewise model 2 
excluded strain constraints associated with another set of elements. Generation of an optimum solution for 
model 1, for example required about 39 CPU minutes. Optimum weight for model 1 was 16 percent 
lighter than that of the original design. The weight savings was 20 percent for model 2. Design changed 
throughout the structure except for variables 12 and 13, which represented the minimum thicknesses. 
Maximum reduction was observed for design variable 8 with a 44 percent change. The thickness reduced 
to 0.44 in. if the original value was 1 in. In other words, in the original structure the entire region that was 
associated with design variable 8 represented an overdesigned condition. The thickness variable 8 for 
model 1 increased to 124 percent from its assumed value of 100 percent, or the region originally was 
underdesigned. The design process redistributed the strain field with many active constraints. There was 
little change in the displacement state or frequency value for model 1. For model 2, maximum 
displacement reduced by 2.5 percent, but its frequency increased by 25 percent because of a 20 percent 
reduction in its weight. 
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Probabilistic Analysis 


Probabilistic analysis requires a geometrical model, a load model, and a material model: 

Geometrical model: All 13 design variables of the deterministic design optimization were considered 
to be random variables. Their mean values were equal to the deterministic optimum design solutions. 
Their standard deviations were set at 7.5 percent of the mean values. The thicknesses of the honeycomb 
were considered as random variables. Their mean values were equal to that of the initial design with a 
standard deviation of 10 percent of the mean values. Cross-sectional areas of the bars were considered to 
be deterministic parameters and were equal to the optimum solutions. 

Load model: Each load component is considered to be a random variable with a mean value equal to 
66.67 percent of the deterministic value, with a 10-percent standard deviation. For example, if a 
deterministic load component is 100 lbf, then its random counterpart would have 66.67 lbf for its mean 
value with 6.67 lbf for the standard deviation. Each load component was changed in a similar manner. 
This refers to load model A. A second load model was recommended and added. Thus, there were two 
load models, load model A and load model B. Load model B was generated following the procedure for 
load model A but with a reduced mean value of 0.50 percent of the deterministic load value. A standard 
deviation was set at 10 percent of the mean load value. 

Material model: Modulus of elasticity as well as Poisson’s ratio were considered to be random 
variables. The standard deviations for all the material parameters were set at 7.5 percent of their mean 
values. The mean values were set to 105 percent of their deterministic values for the Young’s modulus and 
shear modulus of fabric as well as the shear modulus for the honeycomb. The mean value for Poisson’s ratio 
was set to 100 percent of its deterministic value. 


Probabilistic Solution for Strain and Displacement 

The probabilistic solution was obtained using the MSC/Nastran code and the FPI module of the 
NESSUS software. Four different types of distributions were considered: normal distribution, Weibull 
function, lognormal, and Gumbel type 1 distributions. The probabilistic analysis produced the mean value 
of strain (in the element number 11531) to be 5050 microstrain by all four types of distribution functions 
(see row 1 in Table VII). Likewise, the standard deviation remained the same at 16 percent of the mean 
value by all four distribution functions. The mean value of displacement (at node 1 1243 in the 
z-direction) was 9.6 in. with a standard deviation of 10 percent by all four distributions. The observation 
that the mean values and standard deviations did not change for different distribution types was along the 
expected line. 


TABLE VII.— PROBABILISTIC RESPONSE FOR DIFFERENT PROBABILITY LEVELS 3 


N samples 

Micro strain in element 11531 

Displacement at node 1 1243 
along z-direction 


Normal 

Weibull 

Lognormal 

Gumbel 
type 1 

Normal 

Weibull 

Lognormal 

Gumbel 
type 1 

Mean value 

5050 

5050 

5050 

5050 

9.7 

9.7 

9.7 

9.7 

Standard deviation 

810 

810 

810 

810 

0.98 

0.98 

0.98 

0.98 

2 

5050 

5123 

4986 

4917 

9.7 

9.7 

9.6 

9.4 

1000 

7555 

6967 

8162 

9050 

12.6 

11.7 

13.1 

14.4 

50 000 

8379 

7398 

9599 

11523 

13.6 

12.2 

14.5 

17.4 

100 000 

8507 

7459 

9843 

11961 

13.7 

12.2 

14.7 

17.9 

500 000 

8903 

7591 

10402 

12977 

14.1 

12.3 

15.2 

19.1 

1 million 

9264 

7644 

10641 

13415 

14.2 

12.4 

15.5 

19.7 

1 0 million 

9264 

7804 

11424 

14870 

14.7 

12.6 

16.2 

21.4 


a The first and second rows provide the mean values and standard deviations. Remaining data represent mean values. 
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Probabilistic solution for the strain and the displacement for four different types of distributions for 
different failure rates is given in Table VII. For a 50-percent probability of failure, or one failure in two 
samples, the mean value and the calculated strain are equal at 5050 microstrain for normal distribution. 
The same is not true for Weibull, lognormal, or Gumbel type 1 distributions because of a lack of 
symmetry in such functions. Consider next, 1 failure in 1 million samples. The Weibull prediction was 
conservative. It was about 85 percent of the normal distribution for strain as well as displacement. 
Estimates by Gumbel type 1 distribution was on the higher side. For 1 failure in 1 million samples, it 
predicted about 40 to 45 percent higher strain and displacement than that for the normal distribution 
function. 


Reliability-Based Design Optimization 

Probabilistic optimization is performed for the two deterministic design models, models 1 and 2, as 
discussed earlier (see Table VI). The probabilistic design calculation used the information given for 
stochastic analysis along with additional data required to formulate failure modes or the design 
constraints. It was assumed that mean value of allowable strain was 25 percent higher than its 
deterministic limit, while the standard deviation was set at 7.5 percent of the mean value. For example, 
the limits for strain were set with a mean value of 8 microstrain and a standard deviation of s/10. From 
deterministic optimization calculations, it was observed that only the displacement limitation along the 
z-direction and rotation had some influence in the design while other stiffness constraints were quite 
passive. Thus, for probabilistic design only two stiffness constraints were retained. At node 1 1243, along 
the z-direction, the mean value and standard deviation for the displacement limit were 1.25A and 
0.125A in., respectively. The mean value and standard deviation for rotation limit were set at 1.250 and 
0.1250, respectively. For the design calculation normal distributions were assumed for all random 
parameters. 

The optimization calculation required continuous running of the code for more than 5 days, but the 
execution was smooth and eventless. The stochastic optimum solution was quite similar to that obtained for 
deterministic optimization depicted in Table VI. There were numerous active strain, displacement, and 
rotation constraints. Normalized optimum weight obtained for different failure rates is given in Table VIII. 


TABLE VIII.— PROBABILITY OF FAILURE VERSUS WEIGHT 


Samples, 

N 

Normalized optimum weight | 

Design model 1 and 
load case A 
(strength + stiffness) 

Design model 1 and 
load case A 
(strength only) 

Design model 1 and 
load case B 
[strength or 
(strength + stiffness)] 

Design model 2 and 
load case A 
(strength + stiffness) 

2 

64.68 

64.68 

62.24 

63.20 

10 

67.43 

67.43 

63.57 

64.36 

100 

70.47 

70.47 

65.51 

66.97 

1000 

73.75 

73.75 

66.97 

69.57 

10 000 

76.64 

76.64 

68.23 

72.48 

100 000 

79.28 

79.03 

69.70 

76.25 

1 000 000 

92.93 

82.27 

71.57 

91.21 

1 250 000 

94.88 

82.55 

71.76 

93.41 

2 000 000 

100.00 

83.15 

72.16 

98.60 


The normalized optimum weight was set to 100 units for strength and stiffness constraints for load 
case A in design model 1 for 1 failure in 2 million samples. This design exhibited nine active constraints 
consisting of eight strain and one displacement limitation. For an active stochastic constraint the first 
factor in Equation (7) (based on mean values) was equal to the second factor (which was a function of 
standard deviation and the phi critical parameter) in magnitude but with opposite sign. The frequency was 
15.94 Hz for the optimum solution. One SDO run with 61 p levels required 128 h, or 5.33 days. 
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The Inverted-S-Shaped Graph 

The weight versus probability of success given in Table VIII is plotted in Figures 10 and 1 1. In Figure 
10, the x-axis represents N (as in one failure in N samples) and it begins at N= 2, or 50 percent 
probability of success. This graph represents one-half of the inverted-S graph because probability less 
than 50 percent is not included. This graph is for load case A, design model 1; both strength and stiffness 
constraints are considered, and column 1 in Table VIII contains the weight information. 

The same information is replotted in Figure 11, but a logarithmic scale, or log(7V), is used along the 
x-axis. This figure can be approximated into two linear segments. At the intersection point (SD) both 
strength and stiffness constraints are active. From the origin to the point SD, only strength constraints are 
active. From point SD onwards, both strain and stiffness constraints are active. 



Both stress 



Figure 11. — Inverted-S graph in log scale for design model 1 and load case A. 
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The variation of weight shown in Table VIII was as expected. The weight is least for a 50-percent rate of 
success, or for N = 2. The weight increased when failure rate was reduced. The weight in the first two columns 
in Table VIII coincided up to 1 failure in 10 000 samples because only strength constraints were active. The 
weight increased when both strength and stiffness constraints became active (see column 2, in Table VIII). 
Weight shown in column 4 was lighter than that in column 2 because load model B was less severe than load 
model A. Weight shown in column 5 was lighter than that in column 2 because design model 2 was obtained 
by removing few severe strain constraints from model 1 . 

The distribution of the strain state in the wingtip is illustrated in Figure 12 for three cases. The first 
case is the initial design that was obtained by the industry. The second one is the deterministic solution 
while the third one represents the stochastic design. The strain distribution is uneven for the initial design, 
and strain exceeds 4000 microstrain at a few localized locations. The distribution of the strain state is 
more even for the deterministic as well as the stochastic design solutions. 

Optimum weight for six design cases are depicted in Figure 13 in a bar chart. The first case with a 
weight assumed at 100 units represents the original design. All five calculated design cases have lower 
weight than the original design. The deterministic optimum solution (for strength and stiffness constraints 
with design model 1 and load case A) generated a 15. 5 -percent lighter design than the original design. The 
weight decreased further to 80.9 units when the stiffness constraint was relaxed (see design case 3). 
Stochastic optimization produced a mean value for the weight of 83 units with a standard deviation of 
1.3 units for strength constraints for design model 1 and load case A for 1 failure in 1 million samples (see 
design case 4). The mean value for the weight and a standard deviation increased to 94 and 1 .26 units when 
stiffness constraint was added to strength limitation (see design case 5). The mean value of weight and 
standard deviation reduced to 91.7 and 1.23 units, respectively, for strength and stiffness constraints for 
design model 2 and load case B (see design case 6 in Fig. 13). Overall, weight can be reduced up to 
20 percent depending on the choice for model for design, load, and rate of failure. The standard deviation 
was small and remained at about 1 percent for the three design cases 4 to 6. 

Sensitivity analysis was performed for the principal strain (in element 1 1531) for deterministic as 
well as SDO methods and are depicted in Figure 14. The element strain was most sensitive to load and 
elastic modulus E for fabric as well as for design variable 6 (DV6, see Table VI) that contained element 
11531. The other 19 random variables were not sensitive to the strain. Both deterministic as well as 
stochastic methods identified the same set of variables, namely DV6, load, and Young’s modulus. 



Industrial design: 
Weight = 100 units 



Deterministic design 
optimization: 
Weight = 84.5 



Stochastic design 
optimization 
for one failure in 
2 million samples: 
Weight = 83 


Figure 12. — Optimization process more evenly distributed the strain field in wingtip. 
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Weight normalized to 100 units 
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Case 1: Original design 

Case 2: Deterministic optimum design for design model 1 and load case A 
Case 3: Deterministic optimum design for design model 2 and load case B 
Cases 4 and 5: Stochastic design for design model 1 and load case A 
Case 6: Stochastic design for design model 2 and load case B 
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displacement displacement 
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Figure 13. — Optimum weight for six design cases. 
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Figure 14. — Sensitivity analysis for principal strain in element 11531, where DV6 refers 
to design variable 6 and E is elastic modulus, (a) Deterministic design optimization, 
(b) Stochastic design optimization. 
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The CPU time to solution is given in Table IX. The calculation used CometBoards, which is NASA 
in-house software, along with MSC/Nastran version 2005.5.0 (2005R3), and the FPI module of NESSUS 
level 6.2 code (1995). Calculations used a Red Hat Linux 2.6.9-67.ELsmp O/S, with x86_64 architecture, 
2600 MHz, 4 cpu, 8 GB of memory, and 32-bit numeric format. 

One static analysis cycle required about 5 CPU seconds. The run time increased to 51 s for dynamic 
analysis. Stochastic analysis required about 47 min. Deterministic optimization required about 39 min. 
The CPU time for stochastic optimization was enormous at 126 to 128 h of continuous calculations. 


TABLE IX.— CPU TIME FOR DESIGN AND ANALYSIS OF WINGTIP a 


Activity 

CPU time 

One static analysis cycle in MSC/Nastran 

5 s 

Dynamic analysis to calculate 20 frequencies in MSC/Nastran 

51 s 

Deterministic optimization for design model 1 

39 min 

Stochastic analysis 

47 min 

Stochastic optimization: design model 1 and load case A 

128 h 

Stochastic optimization: design model 2 and load case A 

126 h 


a CPU is central processing unit. 


Conclusions 

The testbed CometBoards with MSC/Nastran and a Fast Probability Integrator successfully generated 
reliability-based design optimization for an industrial strength raked wingtip structure of a Boeing 
767-400 extended range airliner made of metallic and composite materials. The optimization run required 
128 h of execution. 

The optimization methodology requires probabilistic models for load, material properties, failure 
theory, and design parameters. Accuracy of the design solution depends on the models. 

Stochastic optimum weight versus reliability traced out an inverted-S-shaped graph. Weight increased 
when risk was reduced and vice versa. The inverted- S graph degenerated into linear segments when a 
logarithmic scale was used for the x-axis. 

The optimization process redistributed the strain field in the structure and achieved up to a 20-percent 
reduction in weight over traditional design. Optimum weight was comparable between a deterministic 
solution and a highly reliable stochastic design. Inclusion of standard deviation as an additional objective 
function should be examined further because it offers very little variation and has insignificant impact. It 
may be a manufacturing issue with little relevance to analytical design calculations. 

Design of an aircraft structural component has to be obtained for multiple load conditions. The 
maximum strain energy criteria identified a few critical load cases from the many load conditions. The 
design generated for the critical load was satisfactory. This approach reduced calculations to a very great 
extent. 

The design sensitivity can identify critical zones for redesign consideration. Both deterministic and 
stochastic concepts identified identical zones. There was no preference to either concept. 

References 

1. Guptill, J.D., et al.: Extension of Optimization Test Bed CometBoards to Probabilistic Design. 6th 
Annual FAA/Air Force/NASA/Navy Workshop on the Application of Probabilistic Methods to Gas 
Turbine Engines, Solomons Island, MD, 2003. 

2. Wei, X.: Stochastic Analysis and Optimization of Structures. Ph.D. Dissertation, Univ. of Akron, 
Akron, OH, 2006. 

3. Thacker, B.H., et al.: Probabilistic Engineering Analysis Using the NESSUS Software. Struct. Saf., 
vol. 28, nos. 1-2, 2006, pp. 83-107. 

4. Guptill, J.D., et al.: CometBoards User Manual Release 1.0. NASA TM-4537, 1996. 

5. Cochran, W.G.: Sampling Techniques. John Wiley, New York, NY, 1977. 


NAS A/TM— 2009-2 15501 


27 




6. Kleiber, M.; and Hien, T.D.: The Stochastic Finite Element Method, John Wiley & Sons, New York, 
NY, 1992. 

7. Cambou, B.: Application of First-Order Uncertainty Analysis in the Finite Element Method in Linear 
Elasticity. Proceedings of the 2nd International Conference Appl. of Statistics and Probability in Soil 
and Struct. Engrg., London, England, 1975, pp. 67-87. 

8. Cornell, C.A.: First Order Uncertainty Analysis in Soils Deformation and Stability. Proceedings of 
the 1st Conf. Appl. of Statistics and Probability in Soil and Struct. Engrg., 1971, pp. 129-144. 

9. Hisada, T.; and Nakagori, S.: Stochastic Finite Element Method Developed for Structural Safety and 
Reliability. Proceedings of the 3rd International Conference Struct. Safety Reliability, Trondheim, 
Norway, 1981, pp. 409-417. 

10. Handa, K.; and Anderson, K.: Application of Finite Element Methods in Structural Safety and 
Reliability. Proceedings of the 4th International Conference on Structural Safety and Reliability, 
Kobe, Japan, 1981, pp. 385-394. 

11. Benaroya, H.; and Rehak, M.: Finite Element Methods in Probabilistic Structural Analysis: A 
Selective Review. App. Mech. Rev., vol. 41, no. 5, 1988. 

12. Elishakoff, I.; and Ren, Y.J.: Some Critical Observations and Attendant New Results in the Finite 
Element Method for Stochastic Problems. Chaos, Solutions & Fractals, vol. 7, no. 4, 1996, 

pp. 597-609. 

13. Schueller, G.I.: Computational Stochastic Mechanics — Recent Advances. Comput. Struct., vol. 79, 
2001, pp. 2225-2234. 

14. Patnaik, S.N.; Berke, L.; and Gallagher, R.H.: Integrated Force Method Versus Displacement Method 
for Finite Element Analysis. Int. Jnl. Computers Structures, vol. 38, 1991, pp. 377-407. 

15. Patnaik, S.N., et al.: Improved Accuracy for Finite Element Structural Analysis Via an Integrated 
Force Method. Int. Jnl. Computers Structures, vol. 46, 1992, pp. 521-542. 

16. Macsyma, Mathematics and System Reference Manual. Fifteenth ed., Macsyma, Inc., 1995. 

17. ANSYS Release 9.0 Documentation. ANSYS, Inc., SAS IP Inc., 2004. 

18. Shinozuka, M.: Monte Carlo Simulation of Structural Dynamics. Int. Jnl. Computers Structures, 
vol. 2, 1972, pp. 855-874. 


NAS A/TM— 2009-2 15501 


28 















